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Abstract 


In this work, we investigate the evolution of a primordial belt of asteroids, 
represented by a large number of massless test particles, under the gravitational 
effect of migrating Jovian planets in the framework of the jumping-Jupiter model. 
We perform several simulations considering test particles distributed in the Main 
Belt, as well as in the Hilda and Trojan groups. The simulations start with 
Jupiter and Saturn locked in the mutual 3:2 mean motion resonance plus 3 
Neptune-mass planets in a compact orbital conhguration. Mutual planetary in¬ 
teractions during migration led one of the Neptunes to be ejected in less than 10 
Myr of evolution, causing Jupiter to jump by about 0.3 an in semi-major axis. 
This introduces a large scale instability in the studied populations of small bod¬ 
ies. After the migration phase, the simulations are extended over 4 Gyr, and we 
compare the hnal orbital structure of the simulated test particles to the current 
Main Belt of asteroids with absolute magnitude H < 9.7. The results indicate 
that, in order to reproduce the present Main Belt, the primordial belt should 
have had a distribution peaked at ~ 10° in inclination and at ~ 0.1 in eccentric¬ 
ity. We discuss the implications of this for the Grand Tack model. The results 
also indicate that neither primordial Hildas, nor Trojans, survive the instability, 
conhrming the idea that such populations must have been implanted from other 
sources. In particular, we address the possibility of implantation of Hildas and 
Trojans from the Main Belt population, but hnd that this contribution should 
be minor. 

Subject headings: minor planets, asteroids: general - planets and satellites: dynamical 
evolution and stability 



3 


1. Introduction 


The jumping Jupiter model is a model of planetesimal driven migration of the major 
planets in which an instability phase arises from the occurrence of close encounters between 
the planets. During this phase of close encounters, a Neptune-size planet is scattered by 
Jupiter, and may be ejected from the Solar System. This causes Jupiter’s semi-major 
axis to “jump” inwards by a few tenth of au. The model was first proposed to solve some 
limitations of other planetary migration models, but it has not yet been fully tested against 
all the possible constraints imposed by many specific characteristics of the Solar System. 
The purpose of this paper is to investigate the effect of the jumping Jupiter evolution on 
the dynamics of a primordial belt of asteroids. 


Planetesimal driven migration was first proposed by iFernandez fc Ipl (1198411 as the 
outcome from the gravitational scattering of a remnant disk of planetesimals by the 
recently formed major planets. The planets were initially assumed to be in a more compact 
orbital configuration. The exchange of angular momentum between the planets and the 
planetesimals induces a smooth change of the planets’ semi-major axes resulting in a 


properties of the Solar System (e.g. 


Malhotra 

1993; 

Liou & Malhotra 

1997) 


expla in only very few 


19971 1. Alternatively, 


instability migration models have been proposed, where the planetary system evolves 
through a short phase of strong dynamical instability. In general, this evolution does not 
destabilize the major planets (e.g. Jupiter and Saturn do not suffer mutual encounters 
in this model), but it has deep consequences for the evolution of other Solar System 
populations. 


The instabili 


by 


Tsiganis et ah 


V moc 


(1299511 


el was introduced by iThommes et ahl (119991) and later reformulated 
, becoming known as the Nice model. In the Nice model, the 


instability arises from the crossing of mutual mean motion resonances between Jupiter and 
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Satur n. The Nice model was able to expl ain the Late Heavy Bo mbardment f|Gomes et al 
20051) . and the origin of Trojan asteroids flMorbidelli et al.ll2005l) . for example. But it failed 


to exp 

ain other 

properties, like the orbits of the terrestrial planets. This led Morbidelli 

et al. ( 

2009 

) and 

Brasser et al. 

(2009 

) to propose a revision of the Nice model, in which the 


instability arises from the gravitational scattering of Uranus by Jupiter and Saturn. In this 
scenario, the successive close encounters cause sudden changes in Jupiter’s (and Saturn’s) 
semi-major axis, so the model is referred to as the “jumping Jupiter” model. A fundamental 
shortback of the original jumping Jupiter evolution is that , more oft e n tha n not, Uranus or 
Neptune is ejected from the Solar System. More recently. 


Nesvornvi (12011 


and Nesvornv 


& Morbidelli (1201211 attempted to solve this limitation considering systems with five initial 


planets: Jupiter, Saturn, and three ice giants. The fifth hypothetical planet is then ejected. 


leaving a system as we known today. This model has had a grea 


success in explaining the 


origin of the irregular satellites of Jup 

iter ( 

Nesvornv et al. 

2014a). several properties of the 

Jovian satellites 

(Deienno et ah! 

2014: 

Nesvornv et al. 

20141 

d, the origin of Jupiter Troians 

(Nesvornv et al. 

2013 

), and the structure of the Kuiper be! 

t Nesvornv (20151. It has also 


proven to fulhll the constraints from the terrestrial planets and the architecture of the outer 
planets. But it has not yet been tested for the asteroid belt. 

In contrast with many other populations of minor bodies, the main advantage of the 
asteroid belt is that it constitutes a very well known population. Its orbital distribution 
is little affected by observational biases, its size distribution is complete up to very small 
sizes, and its taxonomical distribution is well understood. This provides a number of tight 
constraints on the early evolution of the Solar System. In particular, we know that the 
asteroid belt should have suffered a significant primordial excitation, a significant mass 


depletion, and a signi f icant mixture of taxonom ies. Attempts to exp 


first due to 


Wetherilll ( 19921 ). 


review, we refer the reader to 


Petit et ah! (20(11), an d 


Morbidelli et al 


O’Brien et al 


am 


hese issues were 


( 20071 1. For a detailed 


20151 ). Two results are of particular relevance 
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to the present work. 


First, Walsh et ah fl2011 


2 OI 2 I) proposed that the mass depletion and the primordial 


mixing of taxonomies in the asteroid belt could be the consequence of the ear 


of Jupiter and Saturn in a gas disk (iMasset fc Snell grove 


2001 


Pierens et ah 


V mi gration 


2014h . This 


evolution, which presumably happened prior to the phase of planetesimal driven migration, 
became known as the Grand Tack. The time elapsed between the end of the Grand Tack 
(i.e. end of the gas driven migration) and the occurrence of the dynamical instability is 
poorly constrained. It could be as short as a few Myr and as long as ~ 500 Myr, the latter 
assuming that the instability was responsible for triggering the Late Heavy Bombardment 
about 4 Gyr ago. Since our results have implications for the Grand Tack model, we will 
come back to this issue later. 


Second, Minton and Malhotra ( 2009 


20 111) addressed the stirring of the asteroids 


eccentricities and inclinations during the phase planetesimal driven migration. They 
found that in a smooth migration scenario, the characteristic time of migration should be 
unrealistically fast in order to let resonances sweeping to account for the orbital excitation 
of the asteroids. This result favors instability models alia jumpin g Jupiter, that c o uld 


effectively provide the fast migration rates required. In particular. 


Morbidelli et ah 


fl2ninh 


applied a jumping Jupiter model with only four planets to the asteroid main belt, and 
found that it is possible to reproduce its current distribution in inclinations provided that 
the original distribution had an upper cutoff around 20°. We will discuss this result later in 
the light of our findings. 


Our goal here is to analyze and discuss the significance of the jumping Jupiter 
migration in the evolution of the asteroid belt. We perform simulations that reproduce the 
current orbital distribution of the belt and allow us to constrain its initial distribution, that 
is, at the beginning of the planetesimal driven migration. The paper is divided as follows: 























6 


in Sect. [2], we describe the methodology used in our simulations. In Sect. [3] we present our 
results, discussing the implications for the asteroid main belt, the Grand Tack model, and 
the Hilda and Trojan populations. Finally, in Sect. |l]we summarize our conclusions. 


2. Simulations 

We have run a series of numerical simulations of the evolution of the asteroid belt 
from the epoch before the triggering of the planetesimal driven migration to the present 
days. The asteroid belt is represented by a large number of test particles perturbed by the 
Jovian planets. It is worth stressing that the terrestrial planets have not been included 
in our simulations. The initial orbital distribution of the test particles is uniform, and 
the resulting final distribution is compared to the presently observed asteroid belt. This 
allows us to remap the initial distribution onto any desired non uniform distribution and to 
determine the one that provides the best fit to the current asteroid belt. We have assumed 
a strategy consisting in dividing the whole evolutionary process into four stages or phases, 
that we describe in the following. 


2.1. Phase 0: Before the instability 


This phase is intended to simulate the evolution of a primordial asteroid belt between 
the end of the gas driven migration (Grand Tack) and the occurrence of the jumping Jupiter 
instability. During this phase, the planetary system is constituted by Jupiter, Saturn and 
three Neptune-size ice giants. The planets do not migrate and stay in a mut ual resonant 


config uration that is the outcome of the previous gas driven migration phase fjPierens et ah 


20141 ). In particular, Jupiter and Saturn are locked in the 3:2 mean motion resonance, with 
Jupiter slightly outside of its present orbit. The initial osculating values of semi-major axis 
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a, eccentricity e, and inclination / are shown in Table [T] This configuration is stable over 
the simulation time span of 500 Myr. 

The test particles are grouped into three populations; the Main Belt, the Hilda 
group and the Trojan swarms. Each population is represented by 10,000 initial conditions 
uniformly (randomly) distributed over the following intervals: 


• 1.5 < a < 5.0 an, perihelion distance q > 1.5 an, aphelion distance Q < 5.0 an, 
I < 30°, for the Main Belt (hereafter MB set); 

• 4.0 < a < 4.3 an, e < 0.4, I < 30°, for the Hilda group (hereafter HG set); 

• 5.2 < a < 5.7 an, e < 0.2, I < 40°, for the Trojan swarms (hereafter TS set). 


The initial angles M,u,Q (mean anomaly, argument of perihelion, and longitude of node) 
are randomly distributed between 0 and 360° in all cases. Figure [T] shows the initial 
distribution of the test particles. The simulations are carried out using the SWIFT_RMVS4 
symplectic integratojll, with a time step of 0.05 yr. The test particles are eliminated if they 
approach within the radius of any planet, or if q gets smaller than 1 an, or if a evolves 
beyond 100 an. 


Table [2] shows the number of surviving test particles at two specific times: 10 Myr and 
500 Myr. We have chosen these times because they represent the state of the asteroid belt 
in two interesting limits: (i) when the instability is triggered shortly after the end of the 
gas driven migration, and (ii) when the instability happens at a time compatible with the 
Late Heavy Bombardment. We note that at least half of the MB population is stable, while 
the HG and TS populations are significantly erodeco. This reinforces the idea that these 


^https://www.boulder.swri.edu/~hal/swift.html 

^It is worth noting that, due to the partial overlap of the MB and HG sets, the actual 





latter groups are not primordial but implanted from other sources (e.g. 


Levison et al. 


2009 


Nesvornv et al 


paper (see 


2013h . An analysis of the reason 


Morbidell 


2002 


Robutel &: Bodossian 


or th is erosion is beyond the scope of this 


200911 


2.2. Phase 1: Jumping Jupiter 


The purpose of this phase is to simulate the epoch of the jumping Jupiter instability. 


During this phase, the f ive Jovian planets migr a 


previously developed by 


e acc ording to three specific evolutions 


Nesvornv fc Morbidellil (1201211 . These authors performed realistic 


simulations of migration, where the planets interact with a massive disk of planetesimals 
initially located beyond the outermost planet. In the simulations, the planets’ positions and 
velocities were stored in a hie at 1 yr intervals over a tota l time span of 10 Myr. We have 


not reproduced the simulations of 


Nesvornv &: Morbidellil fl2012fl : rather we have mimicked 


the migrati on by reading the stor ed positions and interpolating them using an approach 


described in 


Nesvornv et al. 


(1201311 . Figure [2] shows the three evolutions we have considered. 


which we refer to as Cases 1, 2 and 3, respectively. The main differences between the three 
cases are: 


• The Ice 7^2 planet (middle ice giant) is ejected in Case 1, while the Ice (innermost 
ice giant) is ejected in Cases 2 and 3; 

• Case 2 produces a large number of encounters between Jupiter and the ejected planet, 
which makes Jupiter to suffer a lot of short jumps. Cases 1 and 3, on the other hand, 
produce less interactions and make Jupiter to experience a few longer jumps; 

• In Case 1, the ejected planet reaches minimum heliocentric distances of ~ 2 au during 


number of survivors in the Hilda region is about 1.5 times larger than in the HG set alone. 































9 


the instability, which means that it penetrates deeply in the asteroid belt. In Case 2, 
the heliocentric distance of the ejected planet gets to ~ 1.5 an very briefly. In Case 3, 
the heliocentric distances do not get much below ~ 3 an. 


We recall that in all three cases, the planets are initially in the same orbital configuration 
given in Table [T] The different evol utions in each case arise f rom t he different parameters of 


the planetesimal disk considered by 


Nesvornv fc Morbidellil (120121) . The net inwards jump 


of Jupiter is ~ 0.3 an in all cases, and the instability always occurs between 5.5 and 6.5 
Myr after the start of the simulations. 


These three instability cases have been previ ously tested against a numbe r of 


constraints. They satisfy the constraints defined in 


Nesvornv &: Morbidellil (120121) . namely. 


the final orbits of the outer planets are similar to the real orbits. For example, the proper 
mode in Jupiter’s orbit is excited to its present value by planetary encounters. All three 
cases also satisfy the terrestrial plan ets constraint in tha t Jupiter’s orbit discontinuously 


evolves during planetary encounters (iBrasser et ah 


20091) . This is needed to avoid secular 


resonances with the terrestrial planets, which would otherwise lead to a disruption of the 
terrestrial planets system. All three cases are equally good in explaining the capture and 
orbit al distribution of Jupiter Troians, including their high orbital inclinations f Nesvornv 


et al. 


20131 ). Case 2 shows a richer history of planetary encounters, which leads to a large 


perturbation of the G alilean satellite orbit s that is difficult to reconcile with the present 


orbits of these moons (iDeienno et al. 


2014j) . Cases 1 and 3, on the other hand, have fewer 


planetary encounters and satisfy the Galilean satellites constraint. 


The initial conditions of the test particles for Phase 1 have been cloned from the same 
initial conditions used for Phase 0. The reason for this is twofold; (i) we want to preserve 
the initial angular phases between the planets and the test particles, which is especially 
critical for the HG and TS populations, and (ii) we also want to increase the number 
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statistics. Note, however, that we have only cloned the initial conditions of those test 
particles that survived Phase 0 after 10 Myr and 500 Myr, respectively. This produces two 
different sets of initial conditions for Phase 1. We refer to them as the “Early Instability” 
(El) set and the “Late Instability” (LI) set, respectively. 

For the El set, each initial condition of Phase 0 that survived 10 Myr has been randomly 
cloned 15 times within an interval Aa = ±0.001 an, Ae = ±0.001 and AJ = ±0.1° around 
the reference orbit. The remaining orbital elements were set equal to those of the reference 
orbit. After the cloning process, we have ended up with a MB population of ~ 102 000 
test particles, and a HG±TS population of ~ 15 000 test particles. For the LI set, we have 
applied a similar method, but only to the MB population. We have not taken into account 
the HG and TS populations due to the small amount of survivors at 500 Myr. We have 
created 21 clones of each initial condition of Phase 0 that survived after 500 Myr, ending 
up with a MB population of ~ 99 000 test particles. Figure [3] shows the initial distribution 
of the FI and LI sets. The main difference between the two sets is the larger depletion at 
a > 3.4 an observed in the LI set. 


The simulations have been carried out using a modihed version of the SWIFT_RMVS3 
symplectic integr ator, that interpolates the stored planetary positions and propagates the 


test particles (cf. 


Nesvorny et ah 


20131 ). The total time span of this phase is 10 Myr, and 


the integration time step is 0.05 yr. The test particles are discarded if they hit any planet, 
or if g < 1 an or a > 100 an. 


2.3. Phase 2: Residual migration 

This phase is intended to simulate the residual migration of the Jovian planets after 
the jumping Jupiter instability. The residual migration is a smooth process caused by the 
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interaction of the planets with the disk of planetesinials. It drives the planets to reach their 
current orbits, making Jupiter to migrate inwards, while Saturn and the two remnant ice 
giants (Uranus and Neptune) migrate outwards. This migration has been simulated using a 
modihed version of the SWIFT_RMVS4 code, that applies a non conservative acceleration 
to each planet. The non conservative acceleration a has the form: 


a = a exp 



r ■ V 

V — 2r] ——r 


( 1 ) 


where t is the time, r, v are the position and velocity vector of the planet, and a, r, 77 are 
constants specihc to each planet. The first term provides a smooth drift in semi-major axis, 
while the second term produces damping in eccentricity. An additional term of the form 
—2C,Vz has been introduced in the 2 ; component of the acceleration to produce inclination 
damping. The values of a, r, 77 , ( were tuned to get each planet in its present orbit at the 
end of the simulation. The initial conditions for both the planets and the test particles 
have been taken directly from the hnal conditions produced by Phase 1. An example of 
the initial conditions is shown in Fig. 01 The most notable feature in this figure is the gap 


around 2 an that has been opened 


ay the up, secular reso nance; this resonance is located 


more or less at its present location (iKnezevic et ah 


199ll ) already at the end of Phase 1. 


As in previous phases, the test particles are discarded if they hit any planet, or if g < 1 
an or a > 100 an. The total time span of this phase is 100 Myr, and the integration time 
step is 0.05 yr. 


2.4. Phase 3: Long term evolution 

This is the last stage of our simulations. It intends to reproduce the evolution of the 
asteroid belt from the end of the planetesinial driven migration to the present days. The 
planets do not migrate any more, and their initial conditions are those resulting from Phase 
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2. The total time span of this phase is 4 Gyr. The evolution has been simulated using the 
SWIFT_RMVS4 code, with a time step of 0.05 yr. 

The initial conditions for the test particles have been taken from the final conditions 
produced by Phase 2. However, at the end of Phase 2 there are typically between 25 000 
and 50 000 surviving test particles, depending on the migration case and on the El or 
LI set. Since it is not possible to simulate such amount of test particles over 4 Gyr in a 
reasonable GPU time, we have performed a down-sampling of the sets. The down-sampling 
process involves the following steps: 


test particles with q < 1.6 an and a < 2.1 an are not considered at all. This is because 
we are not interested at this time in the Mars crossing population (g < 1.6 an), no r 


in th e particles that survive in the Extended Inner Belt (a < 2.1 an; 


Bottke et ah 


2 OI 2 I) . This procedure typically removes ~ 30% of the test particles from the whole 


sample. In principle, these test particles should have been signihcantly depleted 
during previous phases if the terrestrial planets were included in the simulations. 
Actually, the only remnant of this population that we recognize today is the Hungaria 
group, whose analysis is beyond the scope of the present worlj^. 


• test particles with 3.5 < a < 5.3 an are all considered. This is because this region is 
significantly depleted during the previous phases and we want to keep track of the few 
survivors, especially at the HG and TS population (see Sect. 13.5p . Typically, these 
test particles represent less than 0.3% of the whole sample. 


• the remaining test particles with a < 3.5 an are selected at random, producing 
a down-sampled set with only 1 008 initial orbits. This implies that the set is 


^The Hungaria group is above the absolute magnitude limit imposed to our comparison 
population in Sect. [3l 
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down-sampled by a factor /down between 15 and 35, depending on the migration case 
and on the El or LI set. 

Note that this approach assigns a larger weight to test particles with a > 3.5 an than to 
the rest of the simulated particles. Nevertheless, this does not introduce any appreciable 
bias in our results, because these particles contribute very little to the distribution of the 
down-sampled set (see Sect. ED. 

An example of the down-sampling procedure is shown in Fig. El Looking at the cyan 
dots in this hgure, we can see that the main features of the present main belt have been 
already sculpted at the end of Phase 2. In particular, we can clearly see the gaps opened 
by the vq and uiq secular resonances (indicated by dashed lines), as well as the Kirkwood 
gaps opened at the 3:1, 5:2. 7:3 and 2:1 mean motion resonances with Jupiter (indicated at 
the top of the plot). 

During Phase 3, each initial set of 1 008 test particles has been propagated together 
with the four Jovian planets. Depending on the migration case and the El or LI set, 
between 15% and 25% of the particles become discarded by the same reasons as in the 
previous phases. 


3. Results 


We have compared the hnal state of the test particles at the end of Phase 3 with 
the current distribution of the real asteroid belt. We have taken into account the 574 
asteroids with q > 1.6 an, a < 4.2 an, and absolute magnitude H < 9.7, that corresponds 
to diameters larger than 4 0-70 km depending on the assumed albedo. The magnitude limit 


is the same considered by 


Morbidelli et al. 


fj2010l ). The upper cutoff in a has been adopted 


to exclude the Trojan population from the comparison. These asteroids may be considered 
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‘primordial” in th e sense that they 


the Solar System flFarinella &: Davis 


rave n ot suffered catastrophic collisions over the age of 


1992|). They have not suffered signihcant variations of 
their orbital elements either, especially in t erms of senii-maio r axis, because they are too 


big to be affected by the Yarkovsky effect flBottke et ah 


20061) . Moreover, these asteroids 


constitute an observationally complete sample. 


Figure [6] shows the distribution of the 574 real asteroids together with the distribution 
of the 800 test particles with g > 1.6 an and a < 4.2 an from the El set, that survived 
at the end of Phase 3 in migration Case fl .e a,...oo. 

between the two distributions. A more detailed way to compare the results is to look at 
the cumulative distributions in e and /, and the density distribution in a. This is shown 
in Fig. [71 As stated before, we have verihed that the region between 3.5 and 4.2 an does 
not contribute signihcantly to the overall statistics. In fact, this region accounts for only 
3% of the real population considered here, and between 0.5% and 3% of the test particles 
(depending on the migration case and the El or LI sets). 

In Fig. [3 we note that there is a clear excess of high eccentricity (e > 0.25) and high 
inclination (J > 20°) test particles with respect to the real asteroid belt. In particular, 
we have verihed that the high inclination particles had high inclinations already at the 
beginning of Phase 0. Therefore, this excess could be eliminated, for example, by applying 
a suitable upper cutoff to the initial population. More realistic initial distributions can be 
tested to hud the one that best matches the real asteroid belt, as we show in the following. 


^It is worth noting that, at the end of Phase 3, Jupiter is at a ~ 5.1 an, that is a little 
bit inwards than the actual location (a ~ 5.2 an). This produces a shift of all the resonance 
locations (Kirkwood gaps) in our simulated test particles. To correct this, we have applied 
an appropriate shift in the semi-major axis of the particles to get the Kirkwood gaps in their 
right position. 
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Our first goal has been to match the cumulative distributions in e and /; the differential 
distribution in a is considered later (Sect. 13.4p . The procedure involved to remap the 
uniform initial orbital distribution of the test particles into a non uniform distribution by 
selecting appropriate initial orbits. Then, the final state of these selected orbits was used to 
match the observations. For this purpose, we have tested four possible types of non uniform 
initial density distributions; 


• upper cutoff distribution 


p{x) = 


1 if a: < c 


0 if a: > c 
parametrized by the cutoff c (with x = e, / ); 


Rayleigh distribution 


parametrized by the mode 7 ; 


X 


p{x) = — exp 


T 


X 


2y 


Maxwell distribution 


Pi^) = yf^exp 


X 


2/32 


parametrized by the mode \/2/5; and 


Gaussian distribution 

r)i ' t \ = — __ 

a\plT\ ^ V 2cr2 
parametrized by the mean p and the standard deviation a 


1 f {x-pf 

p[x) = —exp 


( 2 ) 


( 3 ) 


( 4 ) 


( 5 ) 


Assuming a given type of initial distribution for e, and a given type of initial distribution 
for / (not necessarily the same), we have scanned all the possible values of the distribution 
parameters. Then, we have applied a Kolniogorov-Smirnoff (K-S) test to compare the final 
cumulative distributions in e and / obtained for each value of the parameters, with the real 
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asteroid belt. The results are presented in the form of color maps, where the color level 
indicates the value of the K-S statistic Hks- In these maps, the horizontal axis gives the 
value of the parameter for the assumed initial distribution in e. The vertical axis gives the 


value of the parameter for the assumed initial distribut ion in /. T 


le white curve encloses 


the region where Dks has the highest significance level fjPress et al 


1992). We recall that in 


the case of upper cutoff, Rayleigh and Maxwell distributions, there is only one parameter to 
vary. In the case of the Gaussian distribution, we have fixed one of the parameters (either 
fi or a) and let the other vary. 


3.1. The Early Instability set 


Figure [8] shows the results from the K-S test applied to the two dimensional distributions 
in the e, I plane, like the one shown in Fig. [6l. We show examples using the initial upper 
cutoff and Rayleigh distributions. Best fits with Dks <0.1 have been found for migration 
Cases 1 and 3. For migration Case 2, it has not been possible to get a good fit because there 
is always an excess of high eccentricity test particles. This indicates that test particles in 
Case 2 always get too much excited in e, but not necessarily in /. It is interesting to recall 
that Case 2 have also had a bad performance in the simulations of the Galilean satellites by 


Deienno et al. 


fj2014l) . because the orbits of these satellites get too much excited. We can 


see that, considering upper cutoff distributions, the primordial main belt should have been 
quite exited in I (up to 20°) and e (up to 0. 2). The upper limit in inclinations is compatible 


with the findings by 


Morbidelli et al. 


fl201Cf) . Considering Rayleigh distributions, the e and 


/ should have been peaked at ~ 0.1 and ~ 10°, respectively. A similar result is obtained 
considering Maxwell distributions. In any case, the best fit values are consistent with 
the mean values of e, / in the current main belt. This result supports the idea that the 
presently excited values of eccentricities and inclinations should have mostly been acquired 
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before the phase of planetesimal driven migration, for example, in the Grand Tack model. 
The jumping Jupiter instability and the subsequent residual migration of the outer planets 
help to disperse even more the already excited e and /, and to sculpt the main belt into its 
current shape (e.g. Fig. [5]). 


3.2. The Late Instability set 

Figure 0 is similar to Fig. [HI but for the LI set. Only results for Cases 1 and 3 
are presented, since once again Case 2 has been unable to provide a reasonable good fit. 
We have not found any signihcant differences with respect to the El set, except for a 
slightly tighter constraint of the best hts. This figure also shows that, in general. Maxwell 
distributions generate tighter constraints than Rayleigh (also observed in the El set). We 
can conclude that, concerning the distributions in e and /, there is no evidence that favors 
the model of Early instability over the Late one, or vice-verse. This is expected since the 
main difference between the El and LI sets concerns the stronger depletion of the outer belt 
(Fig. [3]), beyond the location of the 2:1 mean motion resonance, that has a small weight in 
the whole e, / distributions. 


3.3. Implications for the Grand Tack 


The Grand Tack is a specihc case of gas driven migration of the Jovian planets. It was 
originally proposed to explain the low mass of Mars and the current location of Jupiter 


beyond the ice line fj Walsh et ah 


20111 ). In this model, Jupiter and Saturn formed beyond 


the ice line and started to migrate inwards due to the torque exerted by the gas disk (type 
II migration). Saturn migrated faster than Jupiter, and when Jupiter had reached ~ 1.5 
an, both planets became captured and locked in a mutual 3:2 mean motion resonance. At 
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this point, the interaction of the resonant configuratio n with the gas torque reve rted the 


migration, and both planets started to drift outwards flMasset &: Snellgrove 


2001 1. By the 


time Jupiter reached ~ 5.5 an, the gas had already dissipated and the type II migration 
stopped. 


Durin g the Granc 


outwards flWalsh et al 


Tac 


c, Jupiter crosses the asteroid belt first inwards and then 


201211 . Therefore, this model has deep consequences for the evolution 


of asteroids. In particular, the Grand Tack is thought to be responsible for; (i) the 
significant mass depletion of the asteroid belt, (ii) the primordial excitation of the asteroids’ 


eccentricities and inclinations, and (iii) the partial mixing of taxo nomic classes (especia 


the S-type and C-type) in the main belt (for a detailed review see 


Morbidelli et al. 


201511 . 


ly 


According to 


Walsh et al. 


(1201211 . the Grand Tack would produce an asteroidal 


population that can be reasonably well approximated by a Gaussian distribution in e and 
a Rayleigh distribution in I. The best fit Gaussian distribution has mean ^ ~ 0.38 and 
standard deviation a ~ 0.17, while the best fit Rayleigh distribution has mode 7 ~ 10°. 
Using these values, we may therefore test whether the initial distributions predicted by the 
jumping Jupiter instability are compatible with the expectation for the distributions from 
the Grand Tack. 


Figure [To] shows the result. This figure is similar to Figs. [8]and[9l except that the axes 
correspond to initial distributions like those of the Grand Tack. This figure corresponds to 
the El set; a quite similar result has been obtained for the LI set. Once again, migration 
Gase 2 (not shown here) has been unable to provide a good fit due to an excess of high 
eccentricity test particles in the final distribution. 

We note that, concerning the inclinations, good fits to the current asteroid belt are 
obtained for a Rayleigh 7 ~ 10° that is compatible with the Grand Tack nominal value. 
For the eccentricities, good fits are obtained with /i ~ 0.1 and a ~ 0.1-0.3. However, if we 
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force /i to have a large value, we need also a large value of a > 0.35. On the other hand, if 
we force a to have a moderate value, we ne ed a small va l ue of u < 0.15. This result is not 


compatible with the nominal values derived 


Walsh et ah 


(120121) . Actually, our simulations 


using the nominal values produce a hnal distribution with a clear excess of high e particles, 
as shown in Fig. [TTl 


There are some possible explanations for this problem with eccentricity: 


• the distribution produced by the Grand Tack is correct, but the population becomes 
signihcantly depleted at the high eccentricities (e > 0.2) by some mechanism not 
accounted for in our simulations. This depletion might be caused, for example, by the 
terrestrial planets; 

• the distribution produced by the Grand Tack is not correct, and either the model 
overestimates the stirring of eccentricities (i.e. p is too large), or underestimates its 
spreading (i.e. a is too small); 

• the Jovian planets evolved in a different way than considered in our study, thus 
leading to different results than presented here. 

An analysis of these possibilities is beyond the scope of this paper, but the possible effect 
of the terrestrial planets is worth of discussion. 

On one hand, the terrestrial planets might be responsible for a pre-instability depletion 
of the high eccentricity asteroids, that would result in an initial distribution peaked at low 
e and with a moderate spreading, as required (see Fig. [101 niiddle column). This would 
put a strong constraint on the timing of the instability, because: (i) the terrestrial planets 
were almost certainly not completely formed at 10 MyJ^ after the beginning of Phase 0, 

^For example, the Moon-forming impact was dated to have happened at ~ 30-50 Myr. 
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and (ii) even if they would be present at the beginning of Phase 0, it is quite unlikely that 
they produce a signihcant depletion of the belt in only 10 Myr of evolution. In any case, a 
pre-instability depletion model seems to favor a situation where the instability happens at 
later times. Actually, simulations by Deienno et al. (2015, in preparation) indicate that an 
initially cool system of terrestrial planets would require ~ 200 Myr to deplete the main belt 
so as to shift the eccentricity peak of the Grand Tack distribution from 0.38 to ~ 0.2 before 
the instability. These authors also found that the current eccentricities and inclinations 
of the asteroid belt are quite compatible with the Grand Tack initial distribution, except 
at the very low values of e and I. Their model accounted for the terrestrial planets, but 
they considered a simplihed model of the jumping Jupiter instability, in which the orbits 
of the major planets are artihcially moved (instantaneously) from their pre-instability 
conhguration to their present conhguration. 

On the other hand, the terrestrial planets may be responsible for a post-instability 
depletion, since they will certainly continue to erode the main belt at the high eccentricities 
after the instability and until the present times. In this model, the instability may happens 
at earlier times, and the main belt may still loose the necessary fraction of high eccentricity 
orbits so as to reach the present distribution. 

In the end, although including the terrestrial planets might help to conciliate our 
results with the Grand Tack model, they would not allow us to decide between an early vs. 
late instability model. 


3.4. The distribution in semi-major axis 

Once we have determined the best £t parameters for the distribution in e and I, we 
may now pay attention to the distribution in a. Analyzing the histograms like the one in 



21 


Fig. [71 we conclude that our simulations have been able to reasonably reproduce the overall 
distribution in a of the current asteroid belt. Migration Case 1 gives slightly better results 
than Cases 2 and 3. We have not found signihcant differences between the El and LI sets. 
We also realized that the distribution in a is nearly independent on the assumed initial 
distributions in e and I. 


A major discrepancy between our simulations and the real asteroid belt occurs in the 
interval 2.8 < a < 3.0 an. Today, this region of the belt shows a much lower density of 
asteroids compared to th e neighboring region s. The reason for this is still unclear (e.g. 


Minton &: Malhotra 


200911 . 


Broz et ah 


(1201311 refer to this region as the “pristine zone”. In 


our simulations, we have not been able to reproduce the low density of the pristine zone. 
Actually, in some simulations the density we have obtained can be up to twice the observed 
one. 


Here, we test whether a primordial depletion of this zone (i.e. before the jumping 
Jupiter instability) could explain the current density. We do not intend to provide a 
dynamical explanation for such hypothetical depletion, but simply to determine if it could 
be a plausible alternative. The depletion has been simulated by remapping the initial 
density distribution in a into a “square band cut” distribution, i.e.: 

I 1 if O < Q-min or (X > Omax , , 

p(a) = (6) 

£ if Q-niin ^ ® ^ ®max 

with 0<£<1(£ = 0 means that all the test particles within the band are removed). The 
hnal density distribution obtained from the remapped initial distribution has been htted 
to the real asteroid belt applying a test. We have considered a band cut distribution 
with a hxed band width Omax ~ Omin = 0.15 an, compatible with the current width of the 
pristine zone. The band center Oc and the cut level e have been taken as free parameters 
of the fit. Figure IT^ shows the results for the El set in the three migration cases. The 
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color scale gives the values of cind we have considered that good fits are obtained for 
< 0.14 (dark blue). In these examples, we have assumed that e and / initially followed 
the Rayleigh distributions with 7 = 0.1 and 10°, respectively, corresponding to the best fits 
for these orbital elements. 

In Case 1, we observe that e < 0.3 (i.e. > 30% of depletion) around 2.85 an (black 
rectangle) appears to provide slightly better fits to the a distribution. However, the results 
are inconclusive and, in principle, good fits could also be obtained without any depletion 
at all (i.e. e = 1). An example is shown in Fig. [131 This is the same simulation shown 
in Fig. 13 but remapping the initial distributions in e and I only. Apart from the very 
good matches in the e and / distributions, as expected, we see that the distribution in a 
shows a better match in the pristine zone, even without having applied any remapping to 
its initial distribution. On the other hand, in Case 2 no good fit can be found in spite 
of the application of a band cut distribution in a. Finally, in Case 3, an initial band cut 
distribution in a produces a much worse fitting than the original distribution. The reason 
for these differences among the three cases may be related to the amount of dispersion 
in a caused in each migration model. Table |3] shows the median, mean and maximum 
dispersion in orbital elements obtained in the different cases. These values are reflecting 
the behavior of the fifth planet in each simulation, as discussed in Sect. 12.21 (Brasil et ah, 
2015, in preparation). It appears that Case 2 causes too much dispersion in a as to blur any 
structure in the initial distribution. On the other hand. Case 3 does not disperse enough 
the semi-major axes, and any sharp structure in the initial distribution will remain in the 
final distribution. For example, applying a large initial depletion around 2.6 au leaves this 
region almost depleted over the whole simulation (because it is not replenished by particles 
dispersed from the neighboring regions), thus producing the large x^ > 0.25 values observed 
in Fig. [I 2 I 
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The values of 6a shown in Table |3] also allow us to infer that the dispersion caused 
by the jumping Jupiter instability would not be enough to produce a significant mixing of 
the taxon omic classes in the main belt^ as suggested bv the observational evidence I Gradie 


& Tedesco 


1982 


DeMeo fc Carry 


2013 


20141) . This supports the idea that the mixing of 


taxonomies is the consequence of a pre-instability mechanism, as for example, the Grand 
Tack model. 


3.5. The Hilda and Trojan populations 

As discussed in Sect. 12.11 any primordial population of Hildas and Trojans that might 
have existed after the Grand Tack, is strongly depleted already during Phase 0. Even in 
the most optimistic model of an Early Instability, the HG is reduced to ~ 10%, and the 
TS to less than 1% (Table [2]). We have also verified that these remnant populations do 
not survive the jumping Jupiter instability. Of the initially 15 000 test particles in the HG 
and TS populations at the beginning of Phase 1 (red dots in Fig. [3]), less than 1% have 
survived the instability but all of them have been scattered out of the 3:2 and 1:1 mean 
motion resonances with Jupiter. We conclude that in an Early Instability hypothesis, the 
survival probability of a post Grand Tack population is < 10“^ for the HG and < 10“® for 
the TS. These probabilities reduce further by two orders of magnitude in a Late Instability 
hypothesis. 

On the other hand, we have realized that after the instability a number of MB test 
particles have became apparently implanted in the region of the HG and TS populations. 
An example is shown in Fig. [HI top row, where the blue dots correspond to the test 
particles that survived in the HG and TS regions at the end of Phase 1. We have also 
verified that most of these implanted test particles had initial orbits at the beginning of 
Phase 1 with a > 3 au (Fig. [HI bottom row). Since Jupiter is at ~ 5.5 au before the 
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instability, this source region is beyond the 5:2 mean motion resonance with Jupiter, i.e. 
it is equivalent to the current outer main belt. Finally, our simulations also indicate that 
less than 7% of the implanted populations have survived the phase of residual migration 
and are still active at the end of Phase 3 (red dots in Fig. IT^ . This allows us to make 
an estimation of the implanting probability, simply dividing the number of surviving test 
particles at the end of Phase 3 by the total number of test particles with a > 3 an at the 
beginning of Phase 1 (it is worth recalling that the implanted test particles have not been 
down-sampled in the transition from Phase 2 to Phase 3). The estimated probability for 
implanting HG particles is ~ 4 x 10“^ in the most optimistic case (El set, migration Case 
3) and < 4 x 10“^ in the most pessimistic case (LI set, migration Case 1). For implanting 
TS particles, the probability is < 

In order to estimate the current fraction of Hilda asteroids that could have been 
implanted from the MB, we have proceeded as follows: 

1 . we have verihed that, in general, our simulations reproduce quite well the distribution 
of asteroids in the interval 3.2 < a < 3.5 an (see for example Figs. [7t, ITTh and [T3bL 
at least in terms of relative density. This region is presently occupied by the group of 
Cybele, and we will refer to it as the Cybele region; 

2 . we have also verihed that the test particles that have survived in the Cybele region 
have originated in the same region as the implanted Hildas, i.e. our simulated Cybeles 
had a > 3 an at the beginning of Phase 1; 

3. if nc is the number of Cybeles and uh is the number of Hildas that have survived in 
our simulations at the end of Phase 3, then the present fraction of possibly implanted 
Hildas can be computed as: 

, _ riH Nc 

/imp 


/down^C 


( 7 ) 
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where /down is the down sampling factor that has been applied to the Cybele region, 
Nc is the cnrrent number of asteroids in the Cybele region, and Nh is the current 
number of Hilda asteroids (both with if < 11, corresponding to diameter greater than 
~ 40 km for albedo ~ 0.05). 


We computed that /mp < 5% in the simulations with the El set, and /imp < 1% in the 
simulations with the LI set. For the Trojans, these fractions are two orders of magnitude 
smaller. This result has implications for the current distribution of taxonomical classes 
among the Hildas. Assuming that the source of the implanted Hildas is dominated by 
C-type asteroids, as is the case of the current outer main belt, we could expect to find 


today 1 or 2 implanted C-type Hildas with if < 11. Actually, 


the Hilda group is do minated by P-type and D-type asteroids (iDahlgren fc Lagerkvist 


up to this magnitude 


imit, 


1995 


Dahlgren et al 


1997I) IT and only one C-type asteroid is known: (334) Chicago. 


The above results support the idea that current Hildas and Trojans are not primordial, 
but were captur e d from ot her Solar System popu lations during the migration instability. 


Morbidelli et al. 


f|2nn5h and 


Nesvornv et al. 


(12013!) have shown that Trojans can be captured 


from the disk of planetesimals initially located beyond the orbit of Neptune. It is po ssible 


that Hilda asteroids could also be captured from this population (ILevison et al. 


20091 ). but. 


at variance with Trojans, a small fraction (< 5%) of them could have also been captured 
from the outer asteroid main belt. 


4. Conclusions 

In this work, we have investigated the evolution of the asteroid belt during the jumping 
Jupiter instability related to the planetesimal driven migration of the Jovian planets. We 

®See also http://ssd.jpl.nasa.gov/sbdb_query.cgi 
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have simulated the dynamical behavior of test particles from the epoch before the instability 
to the present days. We have considered sets of initial conditions distributed uniformly, that 
have been remapped to different non uniform distributions in order to find a best match 
between the final results and the current asteroid belt. Our conclusions can be summarized 
as follows: 


• we have not found any indication that favors a model in which the instability occurs 
early (soon after the end of the gas driven migration), with respect to a model in which 
the instability occurs later (at a time compatible with the Late Heavy Bombardment 
epoch); 

• of the three models of jumping Jupiter instability that we have tested, the one 
identified as Case 2 has been unable to produce any good fit to the present asteroid 
belt, because the final eccentricities get too much excited; 

• in the other two models (Cases 1 and 3), the best fits indicate that, prior to the 
instability, the asteroid belt should have been quite excited in eccentricities (0.1 to 
0.2) and inclinations (10° to 20°). The main effect of the jumping Jupiter instability 
is only to cause a moderate dispersion of the e, I values. This supports the idea that 
most of the current excitation of the main belt occurred before the instability, as it 
would be the case of the Grand Tack model; 


concerning the upper limit for t h e init ial inclination distribution, we have found the 


same results as 


Morbidelli et al. 


( 120101 ); 


our results predict that before the instability, the asteroid belt should have had a 
Rayleigh (or Maxwell) di stribution in i nclina tions peaked at ~ 10°. This is compatible 


with the value found by 
Grand Tack model; 


Walsh et al. 


(120121 1 for the distribution resulting from the 
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on the other hand, our results predict that the eccentricities of the pre-instability 
asteroid belt sho uld have been peak ed at ~ 0.1, which is smaller than the value of 


0.38 assumed by 


Walsh et ah 


(120121) . This discrepancy could be resolved invoking 


an additional mechanism that would deplete the high e orbits during the transition 
between the Grand Tack and the jumping Jupiter instability. In principle, the 
terrestrial planets, which have not been considered in our simulations, could account 
for such depletion; 


• our model is able to reproduce quite well the density distribution in semi-major axis 
of the main belt up to 3.5 an, except in the the interval between 2.8 and 3.1 an, where 
we have only been able to partially reproduce the very small density observed today. 
This could be indicating that such region became even more depleted by some other 
(not yet understood) mechanism after the jumping Jupiter instability; 

• the median and mean dispersion in semi-major axis caused by the instability does not 
seem to be enough to provide the required mixing of taxonomic classes in the main 
belt; 


• the jumping Jupiter instability is able to implant asteroids from the outer main belt 
into the Hilda group with a very low probability, although enough to explain the 
presence of a few big Hildas belonging to the C taxonomic class. 


We wish to thank the anonymous referee for his helpful comments and suggestions. 
This work has been supported by the Brazilian National Research Council (CNPq) through 
fellowship 312292/2013-9 and grant 401905/2013-6 within the Science Without Borders 
Program, and by NASA’s Solar System Workings program. Simulations has made use of: (i) 
the computing facilities of the Laboratory of Astroinformatics (lAG/USP, NAT/Unicsul), 
supported by FAPESP grant 2009/54006-4 and INCT-A, (ii) the Pleiades Supercomputer 





- 28 - 


of NASA’s High-End Computing Capability, and (iii) the cluster of the Department of 
Astronomy of the National Observatory of Rio de Janeiro, acquired through CAPES grant 
23038.007093/2012-13 



-29- 

REFERENCES 

Bottke, W. F., Vokrouhlicky, D., Minton, D., et al. 2012, Nature, 485, 78 

Bottke, Jr., W. F., Vokrouhlicky, D., Rubincam, D. P., & Nesvorny, D. 2006, Ann. Rev. 
Earth Planet. Sci., 34, 157 

Brasser, R., Morbidelli, A., Gomes, R., Tsiganis, K., & Levison, H. F. 2009, A&A, 507, 1053 
Broz, M., Morbidelli, A., Bottke, W. F., et al. 2013, A&A, 551, A117 
Dahlgren, M., & Lagerkvist, C.-I. 1995, A&A, 302, 907 

Dahlgren, M., Lagerkvist, C.-L, Fitzsimmons, A., Williams, I. P., & Gordon, M. 1997, 
A&A, 323, 606 

Deienno, R., Nesvorny, D., Vokrouhlicky, D., & Yokoyama, T. 2014, AJ, 148, 25 

DeMeo, F. E., & Garry, B. 2013, Icarus, 226, 723 

—. 2014, Nature, 505, 629 

Farinella, P., & Davis, D. R. 1992, Icarus, 97, 111 

Fernandez, J. A., & Ip, W.-H. 1984, Icarus, 58, 109 

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 
Gradie, J., & Tedesco, E. 1982, Science, 216, 1405 

Knezevic, Z., Milani, A., Farinella, P., Froeschle, G., & Froeschle, G. 1991, Icarus, 93, 316 
Levison, H. F., Bottke, W. F., Gounelle, M., et al. 2009, Nature, 460, 364 


Lion, J. G., & Malhotra, R. 1997, Science, 275, 375 



Malhotra, R. 1993, Nature, 365, 819 


Masset, F., & Snellgrove, M. 2001, MNRAS, 320, L55 
Minton, D. A., & Malhotra, R. 2009, Nature, 457, 1109 
—. 2011, ApJ, 732, 53 

Morbidelli, A. 2002, Modern celestial mechanics: aspects of solar system dynamics (London; 
Taylor & Francis), ISBN 0415279399 

Morbidelli, A., Brasser, R., Gomes, R., Levison, H. F., & Tsiganis, K. 2010, AJ, 140, 1391 

Morbidelli, A., Brasser, R., Tsiganis, K., Gomes, R., & Levison, H. F. 2009, A&A, 507, 1041 

Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005, Nature, 435, 462 

Morbidelli, A., Walsh, K. J., O’Brien, D. P., Minton, D. A., & Bottke, W. F. 2015, ArXiv 
e-prints, arXiv; 1501.06204 

Nesvorny, D. 2011, ApJ, 742, L22 

—. 2015, ArXiv e-prints, arXiv;1506.06019 

Nesvorny, D., & Morbidelli, A. 2012, AJ, 144, 117 

Nesvorny, D., Vokrouhlicky, D., & Deienno, R. 2014a, ApJ, 784, 22 

Nesvorny, D., Vokrouhlicky, D., Deienno, R., & Walsh, K. J. 2014b, AJ, 148, 52 

Nesvorny, D., Vokrouhlicky, D., & Morbidelli, A. 2013, ApJ, 768, 45 

O’Brien, D. P., Morbidelli, A., & Bottke, W. F. 2007, Icarus, 191, 434 


Petit, J.-M., Morbidelli, A., & Ghambers, J. 2001, Icarus, 153, 338 



Pierens, A., Raymond, S. N., Nesvorny, D., & Morbidelli, A. 2014, ApJ, 795, Lll 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical 
Recipes; The Art of Scientific Computing, Second Edition, Vol. 1 (Cambridge 
University Press) 

Robutel, P., & Bodossian, J. 2009, MNRAS, 399, 69 

Thommes, E. W., Duncan, M. J., & Levison, H. F. 1999, Nature, 402, 635 

Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 

Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, 
Nature, 475, 206 

—. 2012, Meteor. Planet. Sci., 47, 1941 
Wetherill, G. W. 1992, Icarus, 100, 307 


This manuscript was prepared with the AAS DT[^ macros v5.2. 




32 


Fig. 1.— Distribution of the test particles at the beginning of Phase 0. Black dots represent 
the main belt (MB) population. Red dots represent the Hilda group (HG) centered at ~ 4.1 
an, and the Trojan swarms (TS) centered at ~ 5.4 an. Note that the positions of these 
populations are moved outwards with respect to their present locations, since the 3:2 and 
1:1 mean motion resonances with Jupiter are shifted to larger semi-major axes. The big 
black dot represents the initial position of Jupiter. 
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Fig. 2.— The three cases of jumping Jupiter instability from lNesvorny et al.l (120131 1 that are 
studied in this paper. For each planet, the plots show the semi-major axis and the perihelion 
and aphelion distances. Black corresponds to Jupiter, red to Saturn, and green, blue and 
cyan to the three icy giants. 
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Fig. 3.— Initial conditions for Phase 1: Early Instability set (left) and Late Instability set 
(right). Black dots correspond to the MB population. Red dots correspond to the HG and 


TS populations (only in the Early Instability set). The big black dot represents Jupiter 
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Fig. 4.— The distribution of the MB population in the Early Instability set, at the end 
of Phase 1 (migration Case 1). About ~ 46 000 test particles (of the originally 102 000) 
survived the instability. Their orbits are used as the initial conditions for Phase 2. The big 
dot represents Jupiter. The gap at ~ 2 an is due to the uq secular resonance. 
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Fig. 5.— The tiny dots represent the distribution of the MB population in the Early Insta¬ 
bility set, at the end of Phase 2 (migration Case 1). About ~ 32 000 test particles (of the 
initially 46 000) survived the residual migration. Test particles with q < 1.6 an and a < 2.1 
an (tiny black dots) are not taken into account. Of the remaining ~ 22 000 test particles 
(tiny cyan dots), we perform a random down-sampling to get the set of 1 008 particles repre¬ 
sented by the blue dots (note that particles with a > 3.5 an have not been down-sampled). 
These are used as initial conditions for Phase 3. In red, we indicate the location of the main 
mean motion resonances with Jupiter (left), and the location of the main secular resonances 
(right). 
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Fig. 6.— The blue dots represent the final distribution of 800 test particles in the MB 
population (Early Instability set, migration Case 1), at the end of Phase 3. The distribution 
of 574 real asteroids with H < 9.7 (red dots) is shown for comparison. The dotted-dashed line 
in panel (a) corresponds to g = 1.6 an. The dashed lines in panel (b) give the approximate 
location of the vq and viq secular resonances. 




Fig. 7.— Cumulative distributions in e and /, and histogram of the density distribution in 
a, comparing the real asteroids with 77 < 9.7 (full lines) to the simulated test particles at 
the end of Phase 3 (dashed lines and gray histogram). The simulation is the same presented 
in Fig. [6l 
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Fig. 8.— Comparison of the 2-D distribution in the e, / plane between the test particles of 
the El set at the end of Phase 3 and the current asteroid belt. The color scale gives the 
values of the two dimensional Kolniogorov-Smirnoff statistic T*ks- Best fits are obtained 
for Dks <0.1 (dark blue color). The white curve encloses the region where Dks has the 
highest significance levels. The black pixels in the lower left corner of the plots correspond 
to values for which the K-S test could not be applied due to the small amount of data 
{N < 5). The abscissas give the parameter value of the assumed initial distribution in I 
(i.e. at the beginning of Phase 1). The ordinates give parameter value of the assumed initial 
distribution in e. The top row corresponds to initial upper cutoff distributions. The bottom 
row corresponds to initial Rayleigh distributions. Each column corresponds to a different 
migration case. 
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Fig. 9.— Similar as Fig. |H]but for the LI set. The top row corresponds to migration Case 
1. The bottom row corresponds to migration Case 3. The left column shows initial upper 
cutoff distributions; middle column shows initial Rayleigh distributions, and right column 
shows initial Maxwell distributions. 
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Fig. 10.— Similar to Fig. [SI but assuming initial distributions as those assumed for the 
Grand Tack. In the left column, the abscissas give the mode of an initial Rayleigh distribution 
in /, while the ordinates give the standard deviation of an initial Gaussian distribution in 
e, with hxed mean (/i = 0.38). In the middle column, the abscissas give the mode of an 
initial Rayleigh distribution in J, while the ordinates give the mean of an initial Gaussian 
distribution in e, with hxed standard deviation (a = 0.17). Finally, in the right column, the 
abscissas give the mean and the ordinates give the standard deviation of an initial Gaussian 
distribution in e, assuming an initial Rayleigh dist ribution in / with h xed mode (7 = 10°). 
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Fig. 11.— The same example shown in Fig. [7] (El set, migration Case 1), but assuming initial 
distributions remapped into the nominal Grand Tack distributions (/r = 0.38, a = 0.17, and 
7 = 10 °). 
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Fig. 12.— Results of the f^st (color scale) for the semi-major axis distribution of the 
El set in the three migration cases considered here. The horizontal and vertical axes give 
the values of the parameters Oc and e of the band cut initial distribution in a. A cut level 
£ = 0 means that all the test particles within a band centered at Oc have been removed. A 
cut level e = 1 means that no test particle has been removed. Good fits are obtained for 
< 0.14 (note that is always larger than 0.1). Initially, e and / were assumed to follow 
the Rayleigh distributions with modes 7 = 0.1 and 10°, respectively. The black rectangle 
shown in migration Case 1 indicates a region where the fits are the best. Similar results were 
obtained for the LI set. 
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Fig. 13.— The same example shown in Fig. [7] (El set, migration Case 1), but now assuming 
initial distributions remapped into the Rayleigh distributions in e and I. The modes in e 
and / are 7 = 0.1 and 10°, respectively. No remapping of the distribution in a has been 
applied. Note the larger depletion in the “pristine zone” (2.8 to 3.1 an) compared to Fig. [71 
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Fig. 14.— Top row: Distribution of the MB test particles that became implanted in the 
Hilda region (blue dots) at the end of Phase 1. Gray dots represent the MB population. 
Red dots are the implanted test particles that survived until the end of Phase 3. Bottom 
row: Initial distribution, at the beginning of Phase 1, of the MB test particles that become 
implanted in the Hilda region. Gray dots represent the whole MB initial population. This 


example corresponds to the El set, migration Gase 3. 
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Table 1: Initial orbital elements of the five major planets. 


Planet 

Mass [Mjup] 

a [an] 

e 

/i'i 

Jupiter 

1.0 

5.47 

0.003 

0.05 

Saturn 

0.299 

7.45 

0.011 

0.02 

Ice #1 

0.053 

10.11 

0.017 

0.11 

Ice ^2 

0.053 

16.08 

0.006 

0.07 

Ice #3 

0.053 

22.17 

0.002 

0.05 


Table 2; Statistic of surviving test particles during Phase 0. 


t [Myr] 

MB 

HG 

TS 

10 

6 820 

992 

44 

500 

4 690 

12 

1 


Table 3: The median, mean and maximum dispersion of orbital elements (absolute value of 
the difference between the elements at the beginning of Phase 1 and at the end of Phase 3), 
for each case of migration (El set). 


|5a| [an] \6e\ |5/| [°] 



Median 

Mean 

Max 

Median 

Mean 

Max 

Median 

Mean 

Max 

Case 1 

0.064 

0.114 

0.707 

0.069 

0.112 

0.370 

2.16 

4.25 

22.38 

Case 2 

0.118 

0.231 

1.013 

0.100 

0.137 

0.396 

4.10 

6.56 

26.50 

Case 3 

0.022 

0.072 

0.466 

0.057 

0.095 

0.398 

1.71 

3.16 

19.37 







